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Abstract 

The transportation problem in the plane - how to move a set of objects 
from one set of points to another set of points in the cheapest way - is a 
very old problem going back several hundreds of years. In recent years the 
solution of the problem has found applications in the analysis of digital 
images when searching for similarities and discrepancies between images. 
The main drawback, however, is the long computation time for finding 
the solution. 

In this paper we present some new results by which the time for solving 
the transportation problem in the plane can be reduced substantially. As 
cost-function we choose a distance-function between points in the plane. 
We consider both the case when the distance-function is equal to the 
ordinary Euclidean distance, as well as the case when the distance- function 
is equal to the square of the Euclidean distance. This latter distance- 
function has the advantage that it is integer-valued if the coordinates of 
the points in the plane are integers. 
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1 Introduction. 

The classical transportation problem can be formulated as follows. Let {Sm n — 
1,2, ...,N} denote N sources, let {Rm, rn = 1,2,...,M} denote M sinks, let 
On denote the amount of goods that is available at source Sm let bm denote 
the amount of goods required at the sink Rm, and let c{n,m) denote the 
cost to send one unit of goods from the source Sn to the sink R„i- We call 
a„ the storage at source Sn and bm the demand at sink Rm- Assume that 
J2i On = J2i^ and let A denote the set of iV x M matrices {y{n, m), n — 
1,2, ...,N, TO = 1, 2, ...,M}, such that Y.m=i vin.rn) = a„, Y.n=i = ^™ 

and y{n, to) > 0, n = 1, 2, TV, m = 1,2,..., M. We call a matrix in A a com- 
plete transportation plan. 
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Problem: Find a complete transportation plan X = {x{n,m)} £ A such that 



This problem, the (balanced) transportation problem, is a basic example 
within the theory of optimization theory, and in standard text books one Tisnally 
presents two sohition methods to the problem, namely the simplex method and 
the primal-dual algorithm. The simplex method was developed by G. Danzig 
in the late 1940s. An early presentation of the primal-dual algorithm for the 
special version of the transportation problem called the assignment problem, 
which occurs if TV = M and a„ = 1, for n = 1,2, ...,7V and b-m = 1, for 
m = 1, 2, Af , was given by H. Kuhn in [7]. In this paper Kuhn writes "One 
interesting aspect of the algorithm is the fact that it is latent in work of D. 
Konig and E. Egervary that predates the birth of linear programming by more 
than 15 years (hence the name, the ^Hungarian algorithm^y\ 

In this paper we shall consider the sources and the sinks as points in the 
plane and write {Sn = («„,j„),n = 1,2,...,A^} for the sources and {Rm = 
{Xm,ym),'m = 1,2, ...,M} for the sinks. We shall also assume that the cost- 
function is a distance-function 6{{i,j), {x,y)) between points in the plane. We 
call this special version of the transportation problem for the transportation 
problem in the plane. (From now on we write 6{i,j, x, y) instead of d{{i, j), {x, y)) 
for a distance-function between points in the plane.) 

The transportation problem has a long history, and goes as far back as to 
Monge, 1781 [8]. 

In a paper from 1942 L. Kantorovich proves that the solution of the trans- 
portation problem, can be obtained by solving a maximization problem instead 
of a minimization problem - the so called dual version of the transportation 
problem. (See [6].) He also showed that the solution to the transportation 
problem can be used as a distance-measure between probabilities. 

There is much literature on the transportation problem. Here we only men- 
tion the paper [10] and the book [11] by Rachev in both of which there are many 
further references. 

In the 1980s the solution to the transportation problem in the plane was 
introduced as a distance-measure for digital, grey- valued images with equal total 
grey- value. See [12] and [13]. We call this distance- measure the Kantorovich 
distance for images. 

A drawback with the Kantorovich distance for images is, that in case we 
deal with ordinary grey-valued images, the size of the transportation problem 
becomes quite large which implies that the computation time becomes long. 

Standard algorithms for finding the solution to the transportation problem 
- namely the simplex method and the primal dual algorithm, both have a com- 
putational complexity of order 0{N^) in case the data is integer valued, a fact 
which was pointed out already by Werman et al. in [12]. 

In 1995, Atkinson and Vaidya [2] presented an algorithm for the transporta- 
tion problem in the plane, which, for integer valued problems, has a compu- 
tational complexity of order 0{N^ x log{N) x log{C)), where C denotes the 
maximum of the storages and the demands, in case the distance- function is the 

— metric, and of order 0{N'^'^ x {log{N))^ x log{C)) in case the distance- 
function is the Euclidean metric. However, they did not apply their algorithm 
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to digital images, nor did thoy actually present any experimental results, so it 
is difScult to make comparisons with other methods. 

In a recent paper [5] from 1998 (see also [4]), we described an algorithm for 
computing the Kantorovich distance for images, in case the distance- function 
between points in the plane is chosen to be the — metric or the square of 
the Euclidean metric; the algorithm was based on the primal-dual algorithm 
for the balanced transportation problem, and practical experiments indicated 
a computational complexity of order approximately 0{N'^) for square images 
with N pixels. 

The reason we managed to obtain an algorithm, which in comparison with 
the standard primal-dual algorithm for the balanced transportation problem, 
has a lower computational complexity, was because we found a method by which 
we could reduce the search for the so called admissible arcs. 

However, it was only in the case when the distance-function is the L"^ —metric 
we were able to prove that our algorithm computed the Kantorovich distance 
exactly. In case the distance- function 6{xi,yi,X2,y2) is defined by 

5{xi,yi,x2,y2) = {xi - X2f + {yi -y2f (1) 

we were only able to show that our computation leads to the correct result 
by checking an optimality criterion which exists for the primal-dual algorithm. 
(As we shall motivate in Section 6 below, there are several reasons why it is of 
interest to use a distance-function defined by (1).) 

The purpose of this paper is threefold. One purpose of this paper is to prove 
that the stopping criterion introduced in [5] does indeed lead to the correct 
result in case the distance- function is defined by (1). 

A second purpose is to prove some results by which one can improve the al- 
gorithm for computing the Kantorovich distance for images in case the distance- 
function again is defined by (1). These results imply that we can restrict the 
search for admissible arcs even further. 

One case we did not consider in the paper [5] was the case when the distance- 
function is exactly the Euclidean distance. (We mentioned very briefly, [5] 
Section 20, that in case one uses a linear combination, of the L°° — metric and 
the —metric then the algorithm we described in [5] did work for the examples 
we had considered.) 

The third purpose of this paper is to prove a result by which it is possible to 
improve the primal-dual algorithm for computing the Kantorovich distance also 
in case the distancc-hmction is the Euclidean distance. We believe, that by using 
this result, it should be possible to obtain an algorithm for the transportation 
problem in the plane with Euclidean distance-function, which, for integer valued 
problems, has a computational complexity of approximately order 0{N^) to be 
compared to 0{N^-^) obtained approximately in [2]. 

The plan of the paper is as follows. In Section 2 we introduce some basic 
notations, some terminology and a more precise formulation of the transporta- 
tion problem in the plane. In Section 3 we present the dual formulation of 
the transportation problem and in Section 4 we give a very brief description 
of the primal-dual algorithm. In Section 5 we recall some results proven in 
[5]. In Section 6 we prove the assertion which we presented in [5] in case the 
distance- function is defined by (1), and in Section 7 we prove some other results 
which can be used to decrease the search for admissible arcs when the distance- 
function is defined by (1). In Section 8 we prepare for our results, when the 
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distance- function is the Euclidean distance, by introducing three notions namely 
a hyperbolic set, a level set and an exclusion set, and we prove a simple but very 
important relation between hyperbolic sets and level sets. In Section 9, we prove 
a simple, elementary, lemma which relates hyperbolas and cones, and then in 
Section 10, we apply this result when proving results, by which one can reduce 
the search for admissible arcs in case the distance-function is the ordinary Eu- 
clidean distance. In Section 11, for the sake of completeness, we prove that one 
of the basic assumptions we make in most of our results is correct, namely that 
if one uses the primal-dual algorithm with proper initialization then each pixel 
will always belong to at least one admissible arc. In section 12 we present a 
result by which one can speed up the computation of the quantities by which 
one changes the dual variables when performing the primal-dual algorithm. 

In order for the results of this paper to be useful it is necessary that the 
sources and the sinks can be organized in such a way that one quickly can 
determine all points which are so to speak northeast (northwest, southeast, 
southwest) of an arbitrary point. In section 13 we briefly describe an algorithm 
by which one can accomplish this. Section 14, finally, contains a short summary. 

2 Basic notations and terminology. 

Let K = {{in,jn), n = 1,2, ...,N} be a set of points in R^. By an image P 
with support K (defined on K) we simply mean a set P = {{in, jn,p{in,jn)), n = 
1, 2, .... A^} for which p{in,jn) > 0, n = 1, 2, N. We usually use the notation 
P = {p{i,j) ■ {i,j) G K} for an image. We call an element {i,j) of the support 
K of an image a pixel. 

Next, let P = {p{i,j) : {i,j) E Ki} and Q = {q{x,y) : (x,y) S K2} be 
two given images, defined on the two sets Ki = {{in,jn), n = 1,2, ...,N} and 
K2 = {{xm,ym), m = 1,2, ...,M} respectively. Ki and K2 may be the same, 
overlap or be disjoint. We also assume that 

= X]5(^'y)• 
Let r(P, Q) denote the set of all non-negative mappings h{i,j,x,y) from Ki x 
K2 — > such that 

M 

J2H'i{n)J{n),x{m),y{m)) <p{i{n),j{n)), n=l,2,...,N (2) 

m=l 

and 

N 

^h{i{n),j{n),x{m),y{m))<q{x{m),y{m)), m = l,2,...,M. (3) 

n=l 

We call any function in T{P, Q) a transportation plan from P to Q. A trans- 
portation plan for which we have equality in both (2) and (3) will be called a 
complete transportation plan and we denote the set of all complete transporta- 
tion plans by A(P, Q). 

Let 5{i, i,x,y) be a distance-function measuring the distance from a pixel 
{i,j) in Ki to a pixel {x,y) in K2. The Kantorovich distance dK,s{P,Q) with 
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underlying distance- function d{i,j,x,y) is defined by 

dK,s{P,Q) =min{^ h{i,j,x,y) X S{i,j,x,y) : ■,-,■)£ A(P, Q)}. 

From the definition and some consideration wc sec that computing the Kan- 
torovich distance for images is equivalent to solving a linear programming prob- 
lem called the balanced transportation problem. Most standard text books in 
optimization theory presents an algorithm for solving the general balanced trans- 
portation problem with arbitrary cost-function. See e.g. [1], [3] or [9]. 

3 The dual formulation. 

It is well-known that the solution to the minimization problem described above is 
also obtained by solving the following maximization problem - the dual problem: 

maximize ^ a{i,j)xp{i,j)+ ^ I3{x,y) x q{x,y) 

when 

6{i,j,x,y)-a{i,j)-l3{x,y)>0, {i,j)eKi,{x,y)eK2. (4) 

The variables a(i,j) and P{x,y) are called the dual variables and in case they 
satisfy (4) we have a dual feasible solution. A pair {(i, j), {x,y)} where {i,j) € 
Ki and {x, y) G K2 is called an arc. In case 

o:ii,j)+f3{x,y)=5{i,j,x,y) (5) 

we say that {{i , j) , {x , y)} is an admissible arc. Otherwise the arc is called 

unadmissible. 

If a pixel {i,j) € K\ is such that there exists a pixel {x,y) G K2 such 
that (5) holds then we say that has an admissible arc and vice versa. If 

{(z, j), (x,y)} is an admissible arc we also say that {x,y) is an admissible pixel 
with respect to And vice versa. 

4 The primal-dual algorithm. 

The primal-dual algorithm runs roughly as follows. Let {a{i, j),/3{x, y) : («, j) € 
Ki, {x,y) G K2} be a feasible set of dual variables, and let us also assume 
that initially there exists at least one admissible arc for every pixel in Ki and 
K2. (If this is not the case originally, it is easy to see that one can increase 
some or all of the dual variables so that this hypothesis is fulfilled). Wc now 
look for a transportation plan h{i,j,x,y) from P to Q which has the largest 
total mass among all the transportation plans for which h{i,j,x,y) = in 
case {{i , j) . {x , y)} is unadmissible. In case the transportation plan wc find is 
complete then we are ready; otherwise we update our set of dual variables. To 
find the "best" transportation plan on a given set of admissible arcs we use a 
labeling process, and we use the labeling also in order to find the quantity used 
when updating the set of dual variables. Once we have updated the variables, 
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wc determine the new set of admissible arcs and then we again look for the 
"best" transportation plan on this new set of admissible arcs. Etcetera. 

A good reference on the primal-dual algorithm for the transportation prob- 
lem is Murty [9], chapter 12. Other useful references are [3] and [1]. 



5 Some preliminary results. 

In [5] we were able to improve the primal-dual algorithm for computing the 
Kantorovich distance for images. The reason for the improvement was that we 
were able to rcdiicc the search for new admissible arcs, by using the structure of 
the underlying distance-functions. In this section we shall present some results 
and proofs which in essence can be found in [5] . 

The following proposition holds for any choice of underlying distance- function. 

Proposition 1 Let the underlying distance-function 5{i,j,x,y) be arbitrary. 

Let (ii, ji) and {12, 32) belong to Ki, let [xi. yi) and (2:2,2/2) belong to K2 and 
suppose that {(^2, ^2), (2:2, 2/2)} is an admissible arc. Then 

h) - a{i2, 31) < S{ii,ji,X2,y2) - S{i2,j2,X2,y2) (6) 

and similarly 

/3{xi,yi) - P{x2,y2) < S{i2,j2,xi,yi) - (5(i2, j2, 2^2, 2/2)- (7) 

Proof. Let us prove (7). We have 

P{xi,yi) - 0ix2,y2) = 

P{xi,yi) - 5{i2,j2,X2,y2) + a{i2,j2) < 
5{i2,h,xi,yi) - a(«2,j2) - <5(z2, J2, 2:2, 2/2) + a(»2, j2) = 

<5(«2,i2,2;i,2;i) - S{i2,j2,X2,y2)- 

The proof of (6) can be done in an analogous way. QED. 

By applying the triangle inequality the following proposition follows imme- 
diately from Proposition 1. We therefore state it without proof. 

Proposition 2 Let the underlying distance-function d{i,j,x,y) be a metric. 
Then, if both {ii,ji) S Ki and (12,^2) G Ki have admissible arcs then 

\aiii,ji) - a{i2,j2)\ < S{ii,ji,i2,j2), (8) 

and similarly, if both (2:1,2/1) £ K2 and (2:2,2/2) G K2, have admissible arcs then 

1/3(2:1,2/1) - P{x2,y2)\ < S{xi,yi,X2,y2)- (9) 

Before we state and prove the next lemma let us introduce some further 
terminology. 
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Let (z, j) he a pixel in the support Ki of the image P and let a{i, j) be a dual 
variable corresponding to the pixel (?', j). Let {x,y) be a pixel in the support of 
the image Q. If the dual variable f3{x, y) is such that 



then we say that [x, y) is low with respect to In case there is little risk for 

misunderstanding wc only say that (x, y) is low. 

Let us also introduce the following notation and terminology regarding the 
positions of two pixels. Thus let (xi,yi) and {x2,y2) be two pixels. If xi < X2 
and yi < y2 then wc say that (2:2,2/2) is northeast (NE) of (xi,yi), and that 
{xi,yi) is southwest (SW) of {x2,y2), and if Xi > X2 and yi < y2 then we say 
that {x2,y2) is northwest (NW) of {x\,y\), and that {x\,y\) is southeast (SE) 
of (a;2,?/2)- 

The usefulness of our next result is that it helps to limit the number of 
tests needed for finding all new admissible arcs in case we use the Zi-metric as 
underlying distance-function. The result can be found in [5], (Lemma 19.1) but 

we repeat it here as background information. 

Theorem 1 Suppose that the distance- function we are using is defined by 

S{i,j,x,y) ^\i- x\ + \ j - y\. 

Let (i,j) be a pixel in Ki, which has an admissible arc. Now suppose that 
{xi,yi) e K2, that {x\,y\) has an admissible arc , that {x\,yi) is NE of {i,j) 
and that {xi,yi) is low with respect to {i,j). Then, if {x,y) is NE of {xi,yi), 
then {x.y) is also low with respect to {i,j). (See Figure 1 below.) 



Pix,y) < 6{i,j,x,y) - a{i,j) 



Figure 1 
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admissible 



low 



{xo,yi) 



{xi,yi) 




unadmissible 
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Proof. We prove the theorem by contradiction. Thus suppose that there 
exists a pixel {x, y) S located NE of {x\,y\) and such that at that pixel the 
dual variable (3{x, y) is such that 

-5{i,],x,y) + a{i,]) + (i{x,y) ^ 0. (10) 

But since {xi,yi) is low with respect to (i, j), it follows that 

a{i,j) < 5{i,j,xi,yi) - P{xi,yi) 

which together with (10) implies that 

Sii,j,x,y) - f3{x,y) < 5{i,j,xi,yi) - P{xi,yi) 

and hence 

0{x,y) -P{xi,yi) > S{i,j,x,y) - 5{i,j,xi,yi). 
But since {x,y) is NE of {xi,yi) which is NE of it follows that 

5{i-,j,x,y) - 5{i,j,xi,yi) = x - i + y - j - {xi - i + yi - j) = 

x-xi+y-yi = 5{x,y,xi,yi). 

Hence 

P{x,y) -P{xi,yi) > 5{x,y,xi,yi). 

But since both {x,y) and (xi.yi) have admissible arcs this inequality can not 
be fulfilled because of Proposition 2. QED. 

6 Finding admissible arcs when the distance- 
function is the square of the Euclidean dis- 
tance. 

In [5] we also considered the case when the underlying distance-function is de- 
fined as the square of the Euclidean distance. In our computer experiments we 
first tried to apply the assertion of Theorem 1 in order to reduce the search 
for new admissible arcs, but it turned out, that if we did so, then, for some 
examples, we did in fact exclude too many arcs. Therefore we created another 
claim with slightly stronger assumptions, a claim which we formulated as an 
Assertion, ([5] Assertion 20.1). It turned out that by using this assertion in or- 
der to reduce the search for admissible arcs, we did not exclude any admissible 
arcs from our search. However wc were at that time not able to prove the claim, 
but now we have a proof and can therefore formulate our claim as a theorem. 
We shall first introduce yet some further terminology. 

Consider two pixels (xi.yi) and {x2,y2) in ^^2- Suppose both {xi,y{) and 
{x2,y2) are low with respect to the pixel G Ki. If also 

S{i,j,X2,y2) -a{i,j) - P{x2,y2) > 
S{i,j,xi,yi) - a{i,j) - 0{xi,yi) (11) 

then we say that {x2,y2) is lower than (xi,yi). If we have strict inequality in 
(11) we say that {x2,y2) is strictly lower than {xi,yi). 
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Theorem 2 Suppose that the underlying distance-function is defined by 

5{i,j,x,y) = {i-xf + {j-yf. (12) 

Let {a{i,j),p{x,y)} he a feasible set of dual variables such that each pixel {i,j) 
in Ki, and each pixel {x,y) in K2 have an admissible arc. Let {io,jo) pixel 
in Ki and let {xi,yi), {x2,yi) and {x3,yi) be three pixels in K2 on the same 
line y = yi, such that io < xi < X2 < X3. Furthermore assume that both {xi,yi) 
and (x2,yi) are low with respect to {io,jo) and let {x2,yi) be lower than {xi,yi). 

Then, {x3,yi) is also low with respect to (io,jo)- (See Figure 2 below, for a 
graphical illustration.) 



low lower low 

{X2,yi) {X3,yi) 



{iojo) 

Figure 2. 



Remark 1. This theorem is sUghtly sharper than what we formulate in Asser- 
tion 20.1 of [5] ( sec also [4], Condition 18.1), since we only require that {x2,yi) 
shall be lower than {xi,yi) and not necessarily strictly lower. 
Remark 2. This theorem requires that our images have supports on a regu- 
lar grid structure where it is meaningful to speak of pixels being located on the 
same line. The theorem is not particularly useful for the case when the elements 
of the supports of the two images under consideration are placed at random. 
Remark 3. It ought to be clear that by using this theorem one can decrease 
the number of searches necessary to find all admissible arcs. On each line y = yi 
we do not have to check pixels "further away", once we have found two pixels 
(xi,yi), {x2,yi) such that both are low and (x2,yi) is lower than (xi,yi). 
Remark 4. In Section 7 we shall show how one can obtain an efficient "stop- 
ping criterion" also in the " y-direction" . 

Remark 5. There are several reasons for defining the underlying distance- 
function by (5). It is rotationally invariant, it takes integers into integers, it 
also gives rise to a metric (if one takes the square root afterwards, see e.g. [10]), 
and if the two images are pure translations of each other then in the general 
case there is a unique transportation plan which gives rise to the optimum value. 
Moreover, by comparison with the L^ — metric it seems to give rise to a distance 
which is less course. 
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Before we begin our proof of Theorem 2 we shall prove the following aiixiliary 
result. 

Proposition 3 a) Let (ii,ji) and {12,32) be two pixels in Ki, let {xi,yi) and 
{x2,yi) be two pixels in K2 on the same line and such that xi <X2- Let 6{-, ■,-,■) 
be defined by (12) and suppose that {(ii,ii), {x\,y\)} and {(^2,^2), (2:2,2/1)} are 
admissible arcs. Then ii < 12- 

b) If furthermore {(i2,j2), {xi,yi)} is not an admissible arc then i\ < i2- 
(For a graphical illustration, see Figure 3 below.) 

Figure 3 




Proof of Proposition 3. We first prove part a). The following relations 
must hold: 

oi{ii,3i) + P{xi,yi) = 6{ii,ji,xi,yi) 

Q;(«2,J2) +P{x2,yi) = S{i2,j2,X2,yi) 

"(ii, ji) + f3{x2,yi) < 5{ii,ji,X2,yi) 

a(«2, j2) +P{xi,yi) < S{i2,j2,xi,yi). (13) 

By first adding the last two inequalities and then subtracting the first two 
equalities, and finally shifting terms we obtain the inequality 

S{ii,ji,X2,yi) - d{ii,ji,xi,y-i) > S{i2,j2,X2,yi) - S{i2,j2,xi,yi). (14) 

By using the definition of S{i,j,x,y) (sec (12)), we find that the left hand side 
of (14) becomes equal to 

{x2 - ii)^ - (a;i - iif = xl-x\- 2{x2 - xi)ii 

and the right hand side becomes equal to 

{X2 - i2 f - {Xl - i2f = xl-xj- 2{X2 - Xi)i2. 

Hence in order for (14) to be true it is necessary that i2 > ii- Thereby part a) 
is proved. 
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If furthermore wc know that {(«2,i2), (a^i, J/i)} is not an admissible arc, then 
we have strict inequahty in (13), which imphes that we also have strict inequality 
in (14), and therefore, by using the same kind of arguments as above, we find 
that in order for (14) to be true with strict inequality, it is necessary that i2 > ii- 
Thereby part b) is also proved. Q.E.D. 

Corollary 1 Suppose that the underlying distance-function is defined by (12). 
Let {a(i, j), /3(x, y)} be a feasible set of dual variables such that each pixel {i,j) 
in Ki, and each pixel {x,y) in K2 have an admissible arc. Let {io,jo) belong to 
Ki, let {xi,yi), {x2,yi), {x3,yi) be three pixels in K2 on the same line y = yi, 
such that iQ < x\ < X2 < Xz. Furthermore assume that {xi,yi) is admissible 
with respect to (io.jo) and that {x2,yi) is low with respect to {io,jo)- Then 
{x3,yi) is also low with respect to {io,ja). 

Proof of Corollary 1. Assume that {(iojjo), (2^3, yi)} is an admissible arc. 
Let (12,^2) be an admissible pixel with respect to {x2,yi). Since {x2,yi) is low 
with respect to («o,io), it follows by applying part b) of Proposition 3 to the 
pixels («o, Jo), (*2, ^2), ix2,yi) and {x3,yi), that 12 < io- On the other hand, 
by applying part a) of Proposition 3 to the pixels (io, jo), (12,^2), {x2,yi) and 
(xi, 2/1), we conclude that io < 12 and thereby we have reached a contradiction. 
Q.E.D. 

Remark 1. The corollary implies that in case the pixel (xq, yi) & K2 is the first 
pixel in K2 along the line y = yi for which xq > io and (xq, yi) is an admissible 
pixel with respect to («o, jo), then we can stop the search for admissible pixels 
on this line as soon as we find a pixel which is low. If instead it turns out that 
the first pixel {xo,yi) € K2, with xo > io, along the line y = yi, is not an 
admissible pixel with respect («o, jo), then we should use Theorem 2 until we 
find an admissible pixel, after which we can use the corollary. Clearly the use 
of the corollary only makes a minor improvement of the algorithm. 
Remark 2. Just as is the case with Theorem 2 the corollary is only useful for 
images defined on a grid structure. 

Proof of Theorem 2. We are now ready to prove Theorem 2. Assume that 
{x3,yi) is admissible with respect to (io,jo)- Let (u,ji) be an admissible arc 
with respect to {xi,yi) and let {12, j2) be an admissible arc with respect to 
{x2,yi). From a) of Proposition 3 it follows that ii < 22, from b) of Proposition 
3 we obtain 12 < io, and from the hypothesis of the theorem we also know that 
io < Xi < X2 < X3. Thus 

ii < i2 < io < xi < X2 < X3. (15) 
From the hypothesis of the theorem we also have that: 

a{io,jo) + f3{xi,yi) = S{io,jo,xi,yi) - n (16) 
aiioJo) + fi{x2,yi) = S{io,jo,X2,yi) -r2 (17) 

where 

r2 > ri > 0. (18) 
By subtracting (17) from (16) and using the definition of d{i,j,x,y) we obtain 

(3{xi,yi) - P{x2,yi) = r2 - n + 2{x2 - xi)io + xj - x\ , (19) 
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and by subtracting the equality a{i2,j2) + I3{x2,yi) = S{i2,j2,X2,yi) from the 
inequality a{i2,j2) + f3{xi,yi) < S{i2,j2,xi,yi) we obtain 

fi{xi,yi) - P{x2,yi) < 2{x2 - xi)i2 + xl-xl . (20) 

Finally, by using (19) to replace the left hand side of (20), we obtain 

r2-ri + 2{x2 - xi)io < 2{x2 - Xi)i2 

which implies that 

r\ — r2> 2{x2 — xi){io — 12) > 
where the last inequality follows from (15) and hence 

ri > r2 

which violates (18), and thereby we have reached a contradiction. Q.E.D. 



7 A stopping criterion along the vertical axis. 

Lemma 2 and Corollary 1 imply, roughly speaking, that the search time for 
finding new admissible arcs is decreased from the order 0{N^) to the order 
0{N^-^), in case our images are square digital images, where N as before denotes 
the total number of pixels in the two images. (Lemma 2 and Corollary 1 give 
rise to stopping criteria along the horizontal lines.) We shall now show, how 
we can reduce the; search time to roughly 0{N), by giving a stopping criterion 
along the vertical lines also. 



Theorem 3 Suppose that the underlying distance-function is defined by 

5{i,j,x,y) = {i- xf + {3 -yf- 

Let {Q;(z,j),/3(a;,y)} he a feasible set of dual variables such that each pixel {i,j) 
in Ki, and each pixel {x,y) in K2 have an admissible arc. Let (io, jo) belong to 
Ki. Let (xi,yi) be a pixel in K2, NE of{io,jo), which is not an admissible pixel 
with respect to {io,jo)- Let instead (ii, ji) be an admissible pixel with respect 
to {xi,yi), and assume that also {ii,ji) is NE of{io,jo)- Then if {x2,y2) is a 
pixel in K2, NE of {xi,yi), then (2:2, J/2) is lower than {xi,yi). 
(For a graphical illustration, see Figure 4 below.) 
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Proof. Let us assume that {x2,y2) is NE of {xi,yi) and that it is not lower 
than {xi,yi). Then the following relations hold: 



< r2 < ri. 

By first adding the first equality and the inequality, and then subtracting the 
next two equalities, we obtain 

n -r2 < S{io,jo,xi,yi) + 5{ii,ji,X2,y2) - S{io,jo,X2,y2) - 6{ii,ji,xi,yi), 

and shifting terms we find that 

S{io,jo,xi,yi)-d{io,jo,X2,y2) > S{ii,ji,xi,yi)-S{ii,ji,X2,y2)+ri-r2- (21) 

By using the definition of the distance-function 5{i,j, x, y) we find that the left 
hand side of (21) becomes equal to 



oi{io,3o) + (i{x-i,y\) + r-i = 6{io,jo,xi,yi) 



oi{ii,ji) + P{x2,y2) < S{ii,ji,X2,y2) 
(x{io,jo) + P{x2,y2) + r2 = 6{io,jo,X2,y2) 
Oi{ii,ji) + P{xi,yi) = 6{ii,ji,xi,yi) 



and 



{xi - iof - {X2 - iof + {yi - jof - (2/2 - jof 



13 



x\-xl- 2{xi - x2)io + yl-yl- 2(?/i - y2)jo, 
and that the right hand side of (21) becomes equal to 

{xi - iif - {x2 - hf + (yi - hf - (2/2 - hf + n-r2 = 

xl - x\ - 2(a;i - X2)i\ + Vi - - ^Vi - ^2)^1 + - 

By subtracting the right hand side of (21) from the left hand side, we thus 
obtain the inequality 

2{xi - X2)ii + 2(yi - y2)ii - 2(a;i - X2)io - 2(t/i - y2)jo + r-2 - ri > 0, 

and by shifting terms we obtain the inequality 

2{xi - X2){ii - io) + 2{yi - y2){ji - jo) > n - r2- 

But since x\ < X2, h > io, yi < y2 and ji > jo, and since n > r2 by assump- 
tion, this last inequality can not be true and we have reached a contradiction. 
QED. 

Remark. This theorem gives rise to a very efficient stopping criterion since, 
very often, we can choose the pixel {xi,yi) equal to (ii,ji) equal to (io,2/i) 
which implies that we can stop the search for new admissible pixels on or above 
the line y = yi- 



8 Hyperbolic sets, level sets and exclusion sets. 

Before we start to discuss how we can shorten the search time for new admissible 
arcs in case the underlying distance-function is the Euclidean distance, we shall 
introduce some further notions. As usual, let P = {p{i,j) ■ {i,j) & Ki} and 
Q = {q{x,y) : {x,y) G K2} be two images, let S{i,j,x,y) be an arbitrary 
distance-function, and let {a{i,j) : {i,j) € Kx\ and {(3{x,y) : {x.y) e K2} be 
two sets oi feasible dual variables associated to the images P and Q respectively. 

For each {i,j) e K\ we define the exclusion set E[{i,j),K2] simply as all 
pixels in K2 which are not admissible with respect to {i,j) and similarly for 
each [x, y) G K2 we define the exclusion set E[{x, y), Ki] simply as all pixels in 
Ki which are not admissible with respect to {x,y). 

Next let us also introduce a notion which we call a level set as follows. Let 
(*o, jo) G K\ and let r be a number such that r > 0. The set of all pixels {x, y) 
in K2 such that 

Kio,3o,x,y) - a{io,jo) - P{x,y) > r 
will be called a level set and denoted by L[(io, jo), K2,r\. Note that 

L[{i,3),K2,Q]=E[{i,j),K2]. 

Wc shall now proceed by introducing a notion which \xv shall call a hyperbolic 
set. Let {xi,yi) and {x2,y2) be two points in and let r be a real number. 
We define the hyperbolic set H[{xi,yi), {x2, j/2), r] by 

H[{xi, yi), {x2,y2), r] = {{x, y) & : 5{x, y, X2,y2) - 5{x, y, xi,yi) < r}. 



14 



For a graphical illustration, see Figure 5 below. 




{xi,yi) 

Figure 5. 



Remark. The reason we call the set an hyperbolic set is because if we replace 
the inequality sign in the definition by an equality sign, thereby changing the 
defining inequality to an equality, this equality would define an ordinary hyper- 
bola in case the underlying distance-function is the ordinary Euclidean distance. 

We shall next present a theorem which connects hyperbolic sets and level 
sets. 

Theorem 4 Let P = {p{i,j) ■ & Ki} and Q = {q{x,y) : {x,y) G K2} be 

two im,ages, let 5{i,j,x,y) be an arbitrary distance-function, and let {a{i,j) : 
€ Ki} and {/3{x,y) : {x,y) G K2} be two sets of feasible dual variables 
associated to the images P and Q respectively, such that each pixel in K\, 
and each pixel (x, y) in K2 have an admissible arc. 

Let {io,jo) be a pixel in Ki, let {xi,yi) be a pixel in K2 which is not an 
admissible pixel with respect to (io,jo) ; let the number r be defined by 

r = S{io,jo,xi,yi) - Q:(io, jo) - P{xi,yi), 

let be an admissible pixel with respect to the pixel (xi,yi), let A be defined 

by A = <5(zo, jo) a^i) 2/1) — ^{H,ji,xi,yi), and let s be a real number satisfying 
s > 0. Then, 
(a): 

-ff[(«o,jo),(ii,ji),r- A-s] nK2C L[{io,jo),K2,s], 
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(h): 

H[{io,jo), - A] n if2 C E[{io,jo),K2], 

and (c): 

if also 5{-, ■, ■, ■) is a metric, then 

H[{io,jo), {xi,yi),r- S{io,jo,xi,yi) - s] C i?[(io,.7o), {ii,ji),r - A - s]. 

Proof. Let us first prove (c). We prove (c) by contradiction. Thus assume 
that {x,y) does helongto H[{io,jo),{xi,yi),r — 5{io,jo,xi,yi) — s] but does not 
belong to H[{io,jo), {ii,ji),r- A - s]. Then 

6{x,y,xi,yi) - 5{x,y,io,jo) < r - 6{io,jo,xi,yi) -s 

and also 

S{x,y,ii,ji) - 5{x,y,io,jo) >r- 6{io,jo,xi,yi) + 5{ii,ji,xi,yi) - s. 
Hence 

r-s- 6{io,jo,xi,yi)+5{ii,ji,xi,yi) < 5{x,y,ii,ji) - 6{x,y,io,jo) = 

S{x, y, ii,ji) - S{x, y,xi,yi) + 6{x, y, xi,yi) - 6{x, y, iojo) < 
r - 6{io,jo,xi,yi) - s + 5{x,y,ii,ji) - 6{x,y,xi,yi) 
from which follows after cancellations that 

S{ii,ji,xi,yi) < 5{x,y,ii,ji) - 6{x,y,xi,yi). 

Since we have assumed that 6{-, •, •, •) is a metric, this last inequality can not be 

true because of the triangle inequality. 

In order to prove (a) let (x, y) denote a pixel in K2 and define 

r-2 = 5{iQ,jQ,x,y) - a{io,jo) - f3{x,y). 

Now suppose that the pixel (x, y) docs not belong to the level set L[io, jo, K2, s\. 
This implies that the number introduced above, must satisfy < r2 < s. We 
therefore have the following relations: 

ot{io,3o) + P{x,y) + r2 = 6{io,jo,x,y), 
a{'iiji)+f3{xi,yi) = 5{ii,ji,xi,yi), 
a{io,jo) + I3{xi,yi) +r = S{io,jo,xi,yi), 

and 

a{ii,ji) +P{x,y) < 5{ii,ji,x,y). 

By first adding the first two equations, then subtracting the next equation, and 
finally also subtracting the last inequality, we obtain 

r-2 -r > 6{io,jo,x,y) + 6{ii,ji,xi,yi) - 6{io,jo,xi,yi) - 6{ii,ji,x,y). 

Thereafter shifting terms we obtain 

S{ii,ji,x,y) - 5{io,jo,x,y) > d{ii,ji,xi,yi) - S{io,jo,xi,yi) + r-r2, 
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and using the fact that s > r2 it follows that 

Hii,ji,x,y)-6{io,jo,x,y) > 5{ii,ji,xi,yi)-6{io,jo,xi,yi) + r-s = r-A-s 
from which it follows that 

{x,y) ^ H[{io,jo), {ii,ji),r- A-s] 

which proves (a). 

Finally (b) follows from (a) by taking s = 0. Q.E.D. 
Remark. Part (b) of Theorem 4 can perhaps be considered as the main re- 
sult of the paper. From this result we note that the larger the value r = 
S{io,jo,xi,yi) — a{io,jo) — f3{xi,yi) is, and the smaller the difference A = 
d{iQ,jo,xi,yi) — 5{ii,ji,xi,yi) is, the larger part of the plane will be part of 
the exclusion set. The reason, that we have included part (c) of the theorem, is 
that there are occasions, when it can be simpler to check the size of the hyper- 
bolic set i?[(zo, Jo)) {xi,yi), r — S{io,jo, Xi,yi) — s] than it is to check the size of 
the hyperbolic set H[{io,jo), {ii,ji),r- A - s]. 

9 An auxiliary lemma. 

In this section we shall prove an elementary lemma which relates cones and 
hyperbolas in the plane. We shall rely on this lemma when formulating and 
proving the results in the next section. 

Let 5{-, •, •, •) denote the Euclidean metric. 

Lemma 1 Let a and b be real numbers such that a > and b < 2a, and consider 
the following set: 

HYP = {{x, y)&B? ■ 5{x, y, -a, 0) - 6{x, y, a, 0) > b}. 

Define Hq and H'q by 

Ha = {{x,y)&B? ■.x>Q}, 

H'^ = {{x,y)^R^ ■.x>Q}, 
and for c > define CON E[xo,yQ, c] by 

CONE[xo,yQ,c] = {{x,y) & B? : x > xq, \(y - yo) / {x - xo)\ < c} U (a;o,yo)- 

Then 

(a) : ifb<0 then Hq C HYP, 

(b) : ifb = then H'q = HYP, 

(c) : if < b < 2a and c < ^/{4:a^^ - b'^)/b'^ then CONE[x,y,c] C HYP if 
{x,y)€HYP. 

Proof. Suppose that (x.y) e Ho- Then d{x,y,—a,0) — S{x,y,a,0) > 0, and 
hence if 6 < 0, {x,y) G HYP. This proves (a). 

Suppose that {x,y) e Hq. Then 6{x,y,—a,0) — 6{x,y,a,0) > and con- 
versely if d{x, y, —a,0) — 5{x, y, a, 0) > then a; > 0. Hence ifb = 0, Hq = HYP. 
This proves (b). 
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The proof of (c) is more complicated. Let b> 0. Consider the equation 

d{x, y, —a, 0) — S{x, y, a, 0) = b. 

Taking the square of each side we obtain the equahty 

{x - af + + (a; + af +y^- 2^/{{x - a)^ + y^){{x + a)^ + y^) = h". 

Moving the term 6^ to the left hand side and moving the square root to the 
right hand side and making some simplifications we obtain the equation: 

+ 2y^ + 20^ - 6^ = 2^{{x - af + y^){{x + a)2 + y^). 

Then again taking the square of both sides we obtain 

Ax'^ + V + (20^ - b^f + Sx'^y^ + 8x^a^ - Ax%^ + 8y^a^ - AyH"^ = 

4{{x-af + y^){{x + af + y% 
and the right hand side of this equality can be simplified to 

4x^ + 4a^ - Sx^a^ + 4/ + Sx^y^ + Sy'^a'^. 

Eliminating common terms we obtain the equation 

+ b^- 4a%^ + 8x^a^ - 4x%^ - 4yH^ = Aa^ - 8x^a\ 

and then, making further elimination and moving the terms containing x-factors 
and t/-factors to the left hand side and the other terms to the right hand side, 
we obtain the equation: 

16x^a^ - AxH^ - 4y%^ = -b" + 4a^b\ 

Finally, changing signs and dividing all terms by 46^, this equation can be 
written in a more familiar form as 

y^ = {x^-byA)i{4a'/b')-l). 

This last equation determines a double sided hyperbola, and from the equa- 
tion we also conclude that the lines 

y = ^(4a2 - b^)/b^x 

and 

y = -^(4a2-62)/62a; 

are the asymptotic lines of the hyperbola. Since c < ^J{Ac^^'^^W)lW , it is clear 
from well-known properties of the hyperbola, that the cone CONE\x,y^S^ is a 
subset of HYP as soon as (x, y) € HYP. QED. 
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10 Stopping criteria for the Euclidean distance. 



We shall now consider the problem of how to find new admissible arcs when the 
underlying distance- function S{i,j,x,y) is the Euclidean distance, i.e. 

S{i,j, X, y) = Vii-x^ + U-y)^- (22) 
We shall first introduce some further terminology and notations. We define 
NE[i,j] = {{x,y)GR'' :x>i,y>j} 

and 

LE[i, j] = {{x,y)€R'^ :y = j,x>i}. 

(We have used the letters NE and LE as abbreviations of northeast and "line 
east"). 

Now let (io,jo) be a pixel in Ki. What we are interested in, is to find all 
pixels in K2, which arc admissible with respect to {io,jo)- The purpose of this 
section is to present and prove some results by which one can reduce the search 
time for finding the admissible arcs with respect to {io,jo), which also are NE 
of (*OiJo)- By symmetry it is then easy to reformulate the results so that they 
can be applied, when looking for pixels located NW, SE or SW of {io,jo) which 
are admissible with respect to (io, Jo)- 

In case we have digital images defined on a grid structure, the general al- 
gorithmic procedure is essentially as follows. First check all pixels on the line 
LE[iQ,jo], then on the line LE[io, jo + l], then on the line LE[io, jo + 2], etcetera. 

In order to reduce the search time there are (at least) two ways one can 
accomplish this. Firstly, for each line LE[io,jo + k],k = 0,1,2, .., one can find a 
stopping criterion which implies that one need not to check pixels further away 
from («Oiio) on that line. Secondly, one can find a stopping criterion which 
implies that one does not have to check any of the lines LE[io,jo + k] above a 
certain value of k. 

In this section we shall present two theorems. The first gives stopping criteria 

along the lines LE[io,jo + k],k > 0. The second gives stopping criteria along 
the line {{x, y) ■ x = io,y > jo}- We formulate the theorems in such a way that 
they also can be used when looking for a quantity of the form 

mm{d{i,j,x,y) - a{i,j) - P{x,y) : (z, j) G A, {x,y) G B} 

where A C Ki and B C /1'2- 

However before we state and prove our theorems we shall state and prove a 
simple stopping criterion for pixels on the line iE'fio, jo]- 

Proposition 4 Suppose that the underlying distance-function is defined by (22). 
Let {a{i, j),P{x, y)} be a feasible set of dual variables such that each pixel {i,j) 
in Ki, and each pixel {x, y) in K2, have an admissible arc. Let {io,jo) belong to 
Ki and {xi,yi) belong to K2. Suppose further that {xi,yi) is such that xi > io 
and yi = jo, and that {xi,yi) is not an admissible pixel with respect to (io,jo)- 
Then LE[xi,yi] D K2 C E[{io, jo), K2]. (Recall that E[{i,j),K2] denotes the 
exclusion set with respect to {i,j)). 
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Proof. Our proof follows the same line as our previous proofs. Assume that 
{x,y) G LE[xi,yi] fl K2 and assume also that {x,y) is admissible with respect 
to (zo, jo)- We then have the following two relations: 

S{io,jo,x,y) = a{io,jo) +P{x,y) 

and 

5{iQ,jQ,xi,yi) > a{io,jo) + f3{xi,yi). 
By subtracting the inequality from the equality we obtain 

S{'io,jo,x,y) - d{io,jo,xi,yi) < P{x,y) - f3{xi,yi). 

Since y = yi = jo and x > xi it follows that 

S{io,jo,x,y) - S{io,jo,xi,yi) = x - xi 

and consequently it follows that 

x-xi<p{x,y)-P{xi,y). (23) 

But since 6{-, •, •, •) is a metric and since we have assumed that all pixels have ad- 
missible arcs, it follows from Proposition 2 that (3{x, y)~ (3{xi,y) < 5{x, y, Xi,y) 
and since 6{-, ■) is the Euclidean metric it follows that 

6{x,y,xi,y) = x - xi. 

Hence (3{x, y) — (3{xi, y) < x — xi which combined with (23) implies that x — xi < 
X ~ xi, by which we have reached a contradiction. QED. 

This simple proposition gives rise to a stopping criterion when searching for 
admissible arcs on the line y = jo. Our next aim is to present a result which 
will be useful as a stopping criterion when searching for admissible pixels in the 
sets LE[io,jo + k], fc = l,2, ... . 

Theorem 5 Suppose that the underlying distance-function is defined by 

S{i,j,x,y) ^ Vii-^y + U-y?- (24) 

Let {a{i, j), P{x,y)} be a feasible set of dual variables such that each pixel {i,j) 
in K\, and each pixel {x, y) in K2, belong to an admissible arc. Let (io, jo) belong 
to Ki. Let {xi, yi) be a pixel in K2, NE of (io,jo), and different from («o, jo), 
which is not an admissible pixel with respect to («o, jo)- Let instead (ii, ji) be 
an admissible pixel with respect to {xi,yi), and assume also that (ii, ji) is NE 
of (ioi jo) and that 

ii > io. 

Let r be defined by 

r = S{io,jo,xi,yi) - a(io, jo) - /?(a;i, 2/1), (25) 

let a be defined by 

a = 5{iQ,jQ,ii,ji)/2, (26) 
let s be a real number such that 

< s <r, 
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let b be defined by 

b = S{io,jo,xi,yi) - d{ii,ji,xi,yi) - r + s, (27) 

and let us also assume that b > 0. (The case when b < will be considered in 
the next theorem.) 
Then, if 

> in - jo)/(»i - »o), 

then 

LE[xi,yi]nK2 C L[{io,jo),K2,s]. 

Proof. Let us consider the hyperbolic set H[{io,jo), {ii,ji),r — A — s] where r is 
defined by (25) and where A is defined by A = 6{io,jo,x-i,yi) — S{ii,ji,xi,yi). 
Let the set H be defined by 

H = {{x,y) G : 6{x,y,io,jo)-S{x,y,ii,ji) > b} (28) 

where S{-, •, •, •) is defined by (24) and the number b satisfies (27). Prom the 
definition of if[(zo, jo), (iiiji),''^ — A — s] it is clear that we have 

H = H[{io,jo), {ii,ji),r- A - s]. 

Let us also observe that since < s < r it follows that 

6{xi,yi,io,jo) - 5{xi,yi,ii,ji) > 

S{xi,yi,io,jo) - 5{xi,yi,ii,ji) - r + s = b, 

and therefore {xi,yi) belongs to H. 

Since we have assumed that 6 > and that ii > io, the boundary SH of H 
will be a hyperbola, and the vertex of dH will be along the line between the two 
points (io,io) and and closer to («i, ji) than to (io, jo)- The axis will 

be in the A''£-direction from the vertex, since (ii, ji) is assumed to be NE of 
(«o, jo)- 

In order to guarantee that the set LE[xi,yi] belongs to H it is sufficient 
that the cone Co defined by 

Co = {{x,y) eB? ■.x> XX, \{y-yi)/(x-xi)\ < (ji - jo)/(n -«o)} 

is contained in H. This happens if the angle between the ax;is of the hyperbola 
5H and the asymptotes of 5H is such that the tangent of the angle is larger 
than or equal to (ji — jo)/(H — «o)- We call this angle 9. 

Now by using the estimates obtained in Lemma 1 we conclude that if the 
numbers a and b as defined by (26) and (27) satisfy 

V(4o2 - 62)/62 > {^j^ _ j^y^i^ _ i^)^ (29) 

then the angle 6 between the axis and the asymptotics is so large that any cone 

{{x, y)eR^:x> x' , \{y - y')/{x - x')\ < (ji - jo)/(»i - io)} 

is a subset of H if {x',y') belongs to H. Since we have already proved above, 
that {xi,yi) always belongs to H, it follows that the cone Co is a subset of 
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H. Since evidently LE[xi,y-i] is a subset of Co it now follows from part (a) of 
Theorem 4 that 

LE[xi,yi]n K2 C L[{io,jo),K2,s]. 

QED. 

Proposition 4 and Theorem 5 arc useful for giving stopping criteria along 
lines parallel to the line y = jo- Our next theorem gives rise to stopping criteria 
along the line x = io. The basic assumptions of the next theorem are similar to 
those of Theorem 5. 

Theorem 6 Again, suppose that the underlying distance-function is defined by 
S{i,j,x,y) = ^(i-a;)2 + (j-y)2, 

and let as usual {a{i, j),f3{x, y)} be a feasible set of dual variables such that each 
pixel in K^, and each pixel {x,y) in K2, belong to an admissible arc. Let 

(*0)io) belong to K^. Let (xi, j/i) be a pixel in K2, NE of {io,jo), and different 
from {io,jo), which is not an admissible pixel related to {io,jo)- Let instead 
(iijji) be an admissible pixel related to {xi,yi) and assume that also {ii,ji) is 
NE of (io, jo) (ii not necessarily larger than io). Again, let r be defined by 

r = 6{io,jo,xi,yi) - a{io,jo) - P{xi,yi), 

let a be defined by 

a = (5(jo, jo,ii, ji)/2, 
let s be a real number such that < s < r, and let b be defined by 

b = S{io,jo,xi,yi) - 6{ii,ji,xi,yi) - r + s. 

(a) : Suppose b <0 and that ji > jo- Define 

2/2 = (jo + ji)/2 + {ii - io) V2(ji - jo). (30) 

Then, ify>y2, 

NE[io,y]nK2 C L[{io,jo),K2,s]. 
Next suppose that b > 0, that ji > jo and also that ii > io- Suppose also 

ji-jo<ii-io- (31) 

Then 

(b) : if 

V(4a2 - 62)/62 > (zi - »o)/(ji - jo) 

then 

NE[xi,yi]nK2 C L[{io,jo),K2,s]. 

Next, suppose instead of (31) that j\ — jo > ii— io- 

Then 

(c) : if 

v/(4a2 - 62)/62 > (ji - jo)/(ii - io) 

then 

NE[xi,yi]nK2 C L[{io,jo),K2,s], 
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and 

NE[io, y] n C i[(io, jo), K2, s] 
ify^ys, where the point ys is defined by 

2/3 = ji + 2(ii - iofih - io)/((ji - jof - {ii - iof). (32) 

Proof. Let us again consider the hyperbolic set _ff[(io, jo), (*i,ii), '' — A — s] 
where r is defined by (25), and where A is defined by A = S{io,jo,Xi,yi) — 
S{ii,ji,xi,yi). For each b define the set H{b) by 

H{b) = {{x,y) e : d{x,y,io,jo) - S{x,y,ii,ji) > 6} 

where 5{-, ■,-,■) is defined by (24) and the number b satisfies (27). From the 
definition of H[{io,jo), (*i, ji),?" — A — s] it is clear that we have 

f^(6) = i^[(^o,Jo),(^l,Jl),r-A-s]. 

From the definition of H{b) it is also clear that 

61 < 62 ^ H{b2) C H{bi). 

As we showed in the proof of the previous theorem we also know that {xi,yi) 
belongs to H{b). 

To prove (a) let us first consider the case when 6 = 0. The boundary of H{b) 
will in this case be a straight line, L say, defined by the equation 

y - (.70 + ,7i)/2 = -{{ii - io)/(.n - Jo)){x - {io + h)/2). 

If we consider this equation as a function of a; it is clear, since ii > io and ji > jo, 
that the function is non-increasing, and therefore, since {xi,yi) G H{b), it is 
clear that 

NE[xi,yi]cH{b) (33) 

and also that 

NE[io,y]cH{b) (34) 

if y > 2/2 where j/2 is such that (^0,2/2) are the coordinates of the point where 
the line L cuts the line whose equation is x = ig- To determine the value of j/2 
we just have to insert the value x = io into the equation defining L. We then 
find that 

y = (jo + n)/2 + (in - lofniji - jo), 

and hence y^ satisfies (29). That the assertions of (a) hold, now follows from 
(33) and (34), and hence (a) is proved in case b = Q. But that the assertions of 
(a) are true also in case 6 < follows immediately from the fact that 61 < 62 
if (62) C H{bi). Hence part (a) of the theorem is proved. 

From now on we denote H{b) by H. Next assume that b> and that ii > io- 
Then the boundary of H , which we denote by 5H , will again be a hyperbola, 
the vertex of SH will be along the line between two points (io, jo) and (ii, ji) 
and closer to (ii, ji) than to (io, jo), and the axis will be in the A''i^-direction 
from the vertex, since (ii, ji) is assumed to be NE of (io, jo)- (Recall that H 
is defined by (28)). 
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Wc now want to find conditions on the numbers a and b as defined by (26) 
and (27) sucli tliat the hyperbola dH has an eccentricity which is so large that 
the set NE[xi , yi] is contained in H. Since we have assumed that ji —jo < i\ —io, 
the crucial angle is the angle between the y-axis and the line between the points 
(io,io) and (ii,ji) (that is the line along the axis of the hyperbola). Since we 
have assumed that j\ > jo, this angle is less than 7r/2 and the tangent of the 
angle is clearly equal to {ii — io)/(ii ~ Jo)- By Lemma 1 wc conclude that the 
hyperbola SH will have a sufficiently large eccentricity if the numbers a and b 
will satisfy the equality 

> in - to) /in - jo). 

instead of (29). From Theorem 4 it thus follows that 

NE[xi,yi\nK2 C L[{io,jo),K2,s] 

and thus (b) is proved. 

It remains to prove (c). Thus assume that 6 > 0, that ii > io, and that 
ii — jo > ii — io- Again we consider the set H defined by (28), and again we 
want the eccentricity of the hyperbola SH to be so largo that the set NE[xi,yi] 
is a subset of H. This time the crucial angle is the angle between the a;-axis 
and the line through the two points {io,jo) and («i, ji) (that is the line along 
the axis of the hyperbola). Since we have assumed that ii > io, this angle is 
less than 7r/2, and the tangent of the angle is clearly equal to {ji — jo)/(*i ^ *o)- 
Again, by applying Lemma 1, we conclude that the hyperbola dH will have a 
sufficiently large eccentricity if the numbers a and b will satisfy the inequality 
(29). Hence by Theorem 4 it follows that if (29) holds then 

NE[xi,yi]nK2 C L[{io,jo),K2,s] 

and thus the first part of (c) is proved. 

Moreover, if a and b satisfy the inequality (29), then the line, which starts 
at the point {ii,ji) and which makes an angle with the axis of the hyperbola 
SH which is the same as the angle between the axis of the hyperbola and the 
a;-axis, will be contained in H. This line will cut the line a; = io at some point 
1/3. From geometric considerations it is then clear that NE[io, y,] C H ii y > ys 
and hence by Theorem 4 we also have 

NE[io,y]nK2 C L[io,jo,K2,s] 

if y > ya- 

It remains to show that the number ys is determined by the expression (32). 
The equation for the line we are considering can be written 

y-ji = k{x-ii), 

where the number k still has to be determined. If we denote by 9 the angle 
between the line through the two points (io, jo) and (ii, ji) (that is the line 
along the axis of the hyperbola SH) and the line y = io , then 

k = tan{20) = 2tan{9)/{l - tan^{0)). 
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Now since tan{9) = (ji — — io), by inserting x = iq into the equation for 

the line , we find that 

y = ji + 2(ji - io)(io - n)/(i - ((ji - io)/(n - io)f, 

which after sinipHfication is equal to the expression in (32 ). 

Thereby also the second part of (c) is proved and thereby the proof is com- 
pleted. Q.E.D. 

Remark. By using the theorems of this section, it seems likely, that the search 
time for finding all pixels which are admissible with respect to a given pixel in 
most situations will have a computational complexity of order 0(1), and there- 
fore the computational complexity for finding all admissible arcs will be of order 
0{N) roughly, where though the constant in 0{N) may be fairly large. 

11 Each pixel belongs to an admissible arc. 

One of the basic assumptions we have made in most of the results proven above 
is that the set of dual feasible variables {a{i,j), (5{x, y)} is such that each pixel 
(i,j) S Ki and each pixel (.t,i/) G K2 belong to an admissible arc. That is, to 
each G Ki there exists a pixel {x, y) G K2 such that 

"(ii) + P{^',y) = S{i,j,x,y), (35) 

and to each {x,y) G K2 there exists a pixel {i,j) G Ki such that (35) holds. 
We shall now show that this is a valid assumption independently of the choice 
of the underlying distance- function. 

Let us first state the following proposition. 

Proposition 5 Define the set of dual variables {a{i,j),/3{x,y)} as follows: 

a{i,j) = mm{S{i,j,x,y) : {x,y) G K2} (36) 

P{x, y) = min{5(i, j, x, y) - a{i, j) &Ki}. (37) 

Then the set of dual variables is feasible, and, moreover, to each pixel G Ki 
and to each pixel (x, y) G K2 there exists an admissible arc. 

Proof. The proposition is intuitively obvious. A formal proof can read as 
follows. ^Prom equation (37) it is clear that to each {x, y) G K2 there exists a 

pixel {i,i) G Ki such that (35) holds. Moreover, from (37) it also follows that 
the set of dual variables is a feasible set. Now suppose that there exists a pixel 
G Ki such that 

a(«,j) +/3(a;,2/) < S{i,j,x,y) 

for all {x,y) G i^2- Since a{i,j) = 5{i,j,x,y) for some {x,y) G K2, it follows 
that for this choice of pixel, P{x,y) must be negative. But this is not possible 
because of (37) and the fact that 6{i,j,x,y) > a{i,j) because of (36). Q.E.D. 
We shall now prove: 

Lemma 2 Let the set {a{i, j), /3{x,y)} be a set of feasible dual variables ob- 
tained by the primal-dual algorithm when initially the set of dual variables are de- 
fined as in Proposition 5. Then each pixel {i,j) G Ki and each pixel {x,y) G K2 
will have an admissible arc. 
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Remark. In case one is familiar with the primal-dual algorithm, the truth 
of the lemma is intuitively clear. We include a formal proof for the sake of 
completeness. 

Proof. To give a formal proof of the lemma we need to describe the primal-dual 
algorithm partly again. As pointed out in Section 4, one step in the primal- 
dual algorithm is called the labeling procedure, and below we shall describe it 
in more detail. Let us first introduce two further notions. Thus let h{i, j, x,y) 
be a transition plan from the image P = {p{i,j) ■ € Ki} to the image 

Q = {q{x, y) : {x, y) G K2}. (See Section 2 for the definition of a transportation 
plan.) We call a pixel € Ki deficient if the transportation plan h{i,j, x, y) 
is such that 



If instead (i, j) G Ki is such that the left hand side of (38) is equal to the right 
hand side of (38) then we say that («,j) is full. 

We denote the subset of Ki consisting of all deficient pixels by Di and the 
subset consisting of all full pixels by Fi . 

The labeling procedure is as follows. One starts the labeling procedure by 
labeling all deficient elements in Ki (that is the set Di). Then, whenever a 
pixel G Ky is labeled, then all pixels {x,y) G K2 which arc admissible 

with respect to (i, j) and are not yet labeled, are labeled, and whenever a pixel 
{x,y) G K2 is labeled, then all pixels G K\, which are admissible with 
respect to {x,y), which are not yet labeled and for which h{i,j,x,y) > 0, are 
labeled. 

At the end of a labeling procedure, depending on the labeling, one can 

either go to the so called flow change routine or to the so called dual solution 
change routine. If one goes to the flow change routine the transportation plan is 
redefined after which one can redo the labeling. Sooner or later one has reached 
a state when it is not possible to find a better transportation plan for the given 
set of admissible arcs induced by the present set of dual variables. One therefore 
has to redefine the dual variables. 

When the labeling procedure ends, a number of pixels in both Ki and K2 
are labeled. Let Li denote the set of labeled pixels in Ki, let Ui denote the set 
of unlabeled pixels in Ki, let L2 denote the set of labeled pixels in K2, and let 
U2 denote the set of unlabeled pixels in K2. 

When the dual variables arc changed, one uses the following quantity: 

Q = ma.x{S{x,y,i,j) - a{i,j) - P{x,y) : (i,j) eLi,{x,y) G U2} (39) 

which one can prove will be > (unless the transportation plan is complete). 
The way the dual variables are changed is as follows: 

Oinew{i,j) = Oiold{i,j) + ©, V(i, j) G I/l 

anew{i,j) = aoid{i,j), e Ui 

I3new(x,y) = (3oid{x,y) - 6, V{x,y) G L2 

Pnew{x,y) = (3oid{x,y), y{x,y) G U2. 

From the way the new set of dual variables are defined the following conclusion 
follows immediately. 




h{'. 




(38) 



{x,y)eK2 
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Proposition 6 After we have updated the dual variables the only way by which 
a,n arc {{i,j), {x, y)} can be unadmissible when it previously was admissible is if 
ii,j) e Ui and e L2. 

Now let us thus prove that the new set of dual variables gives rise to a set of 
admissible arcs such that to each pixel in Ki and to each pixel in K2 there exists 
an admissible arc. 

Our induction assumption is that the old set of dual variables has this prop- 
erty. Note that initially this is true because of Proposition 5. Now let us first 
prove that to each pixel {i,j) G Ki there exists a pixel in K2 which is an ad- 
missible arc. From the induction hypothesis and Proposition 6 above, the only 
case we have to consider is the case when is unlabeled. But if is 
unlabeled it must be full. Next let G{i,j) be the set of pixels (x.y) G K2 for 
which {{i, j), {x,y)} is an admissible arc under the old set of dual variables. 
If there exists an element {x,y) e G{i,j) which is not labeled then we have 
found an arc which will also be admissible under the new set of dual variables. 
Otherwise all pixels in G{i,j) must be labeled. Since the pixel is full it 

follows that there must exist a pixel {x,y) G G{i,j) for which h{i,j,x,y) > 0. 
But then because of the way the labeling procedure works the pixel (i, j) would 
have been labeled. Hence if there exists a pixel (i, j) G Ki which is not labeled 
there must exist a pixel {x,y) G G{i,j) which is not labeled and consequently 
this arc will be admissible also after the dual variables are updated. 

We have now proved that to every pixel {i,j) G Ki there exists a pixel 
{x,y) G K2 such that the arc {x,y)} is admissible under the new set of 

dual variables. It remains to prove that each pixel (.x, y) G K2 has an admissible 
arc under the new set of dual variables. But this is somewhat easier. By Propo- 
sition 6, a necessary condition for an arc {{i,j),{x,y)} tohe unadmissible when 
it previously was admissible is that (x, y) is labeled. But from the labeling pro- 
cedure, it then must exist an admissible pixel {i,j) with respect to {x,y) which 
also is labeled, and then because of Proposition 6 the arc {x, y)} will be 

admissible also under the new set of dual variables. Thereby the induction step 
is verified and the proof is completed. Q.E.D 

12 Reducing the computation time for deter- 
mining the quantity by which the dual vari- 
ables are changed. 

An important step in the primal-dual algorithm for the transportation problem 
is to determine the quantity O as defined by (39). In general the computational 
complexity for this part of the algorithm is 0{N^), where though 0{N'^), so to 
speak, has a small constant. 

In case one has integer storages and demands and integer- valued cost-function 
then one can always take © equal to 1 and this brings the computation time 
down to zero. This strategy works also quite well in case we are working with 
digital images. 

In case we are dealing with the transportation problem in the plane and 

have the Euclidean distance as cost-function and the positions of the sources 
and sinks are essentially random, then it is necessary to compute the quantity 



27 



O defined by (39). In this section we shall formulate a result by which one can 
reduce the computation time to determine the quantity 8. 

We first prepare for our result. Let as always Ki = {{i{n),j{n) : n = 
1,2,..., AT} and K2 = {(x(m),y(m) : m = 1,2,...,M} be two sets in R"^, let 
P — {pii,j) ■ G Ki} and Q = {q{x,y) : {x,y) G K2} be two images such 
that p{i,j)=J2K li^^y)- Let the underlying distance-function be defined 
by 

S{i,j,x,y) = V(i - ,x)2 + (j - yy. (40) 

Let {a(i,j), (3{x, y)} be a feasible set of dual variables such that each pixel (i, j) 
in K\, and each pixel (x, y) in K2 belong to an admissible arc. Let L c Ki and 
f/ C be such that 

Q = mi7i{S{i,j,x,y) - a{i,j) - I3{x,y) : G L,{x,y) gU} > 0. (41) 

Our aim is now to determine 6 and we do this in principal by computing 

B(ij) = mm{6{i,j,x,y) - a{i,j) - (3{x,y) : {x,y) G U} 

for each (i, j) G L. 

Let the size of L be equal to A^i, let the size of U be Mi let {{i{k),j{k)) : k = 
1,2, ...,7Vi} be a sequential list of the set L, let {{x{k),y{k)) : k = 1,2,..., Mi} 
be a sequential list of U. For n = 1, 2, A''i, and m = 1, 2, Mi, let 

e(n) = min{e(i(fe) : 1 < fc < n}, 

let 

0(i,i)(w) = mm{d {i,j,x{k),y{k)) - a{i,j) - P{x{k),y{k)), 1 <k< m}, 

let 

9(0, to) = e(i(i)j(i))(TO), 

set 

6(0) = 6(0,1), 

and define 

6(n,TO) =min{6(n),6(i(„+i)j(„+i))(TO)}. (42) 

Theorem 7 Let L = {{i{k), j{k)) : k = 1,2, Ni} C Ki and 

U = {{x{k),y{k)) : k = 1,2,..., Mi} C K2 be such that (4I) holds. Let n and 

m be two integers such that < n < Ni — 1 and 1 < Mi — 1 and n + m > 1 . 

Set (iaJo) = («(n + 1), j'(" + 1)), set {xi,yi) = (.x(to + 1), y(m + 1)), let 

be an admissible pixel with respect to {xi,yi) and assume that (ii, ji) is NE of 

{io,jo). Let r be defined by 

r = S{io,jo,xi,yi) - a{io,jo) - P{xi,yi), 

let a be defined by 

and let b be defined by 

b = 6{io,jo,xi,yi) - 6{ii,ji,xi,yi) - r + Q{n,m) 
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where <d{n,m) is defined by (42). 
Then: 

(a) : Suppose b < and that yi > (jo + ji)/2. Then NE[xi,yi] n K2 C 
L[{iaJa),K2,Q{n,m)]. 

(b) : Next suppose that b > and suppose also that ji > jo- Suppose also 

ji - jo <h-io (43) 

and 

v/(4a2 - 62)/62 > _ i^y^j^ _ j^y 

Then 

NE[xi,yi]nK2 C L[{io,jo),K2,e{n,m)]. 

Suppose instead of (43) that j\ — jo > ii — io > 0. 
Then, if 

then 

NE[xuyi] nK2C i[(io, jo), K2, e(n, m)]. 
Remark. Since the proof is very similar to the proof of Theorem 6, we omit it. 

13 Organizing points along the northeast direc- 
tion 

In order to use Theorems 6 and 7 efficiently, it is necessary to have an initial- 
ization process such that one can find elements NE of a given element rapidly. 
Thus, it is necessary that one can organize the elements in the sets Ki and if 2 
and also the set Ki U K2 in the northeast direction, - as well as the northwest, 
the southeast and the southwest directions. In case we are dealing with digital 
images on a rectangular grid then the pixels are organized already from the 
start in such a way that this is no problem. However, in the general case it is 
necessary to organize the elements of the two given images in the NE — SW 
direction, and the NW — SE direction. In this section we shall describe one 
algorithm by which this can be done. We shall only consider the northeast 
direction since the other directions can be handled analogously. 

Thus let K = {{x{n),y{n)),n = l,2,...,iV} be a set of points in the plane. 
What we want to do is to organize the points in K so that we can apply Theo- 
rems 6 and 7 cfBcicntly. What we will create is a structure which is often called 
a directed acyclic graph. 

We start by introducing an extra point (a;(0), y{0)) such that a;(0) < x{n), n = 
1,2,..., N and y{0) < y{n), n = 1,2, ...,N. 

Next let us define z{n) = x{n) + y{n), n = 1,2, ...,N and let us first order 
the elements of K in such a way that 

z{n) < z{n+l),n= 1,2,..., N - 1 

and 

z(n) = z{n + 1) x{n) < x(n + 1). 
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This way to order the elements of K impUes of course that {x{n),y{n)) can not 
be NE of {x{m),y{m)) if n < m. 

To each element {x(n),y{n)), < n < N we shall associate a number 
of related elements which will be called parents and children. These notions 
are defined as follows: Let {x{n),y{n)) and {x{m),y{m)) be two points in K 
with n < m, and suppose that {x{n),y{n)) is SW of {x{m),y{m)). The point 
{x{n),y{n)) will be a parent of {x{m),y{m)) if cither 1) it is the only clement 
in K which is SW of (x(m), y{m)), or 2) if there are other points in K which 
are SW of {x{m),y{m)) then {x{n),y{n)) is not SW of any such point. If 
{x{n),y{n)) is a parent of {x{m),y{m)), then we say that {x{m),y{m)) is a 
child of {x{n),y{n)). 

In case we define parents and children as above, then, if we have found the 
parents of all points in a set K & E? , then we say that we have organized the 
set in the NE — SW direction. 

It is easy to construct examples of sets of size N for which the number of 
parents will be N'^/A. On the other hand it is clear that one does not need 
more than A''(A'' + l)/2 checks to find all parents. Therefore it is clear that 
the computational complexity to organize a set in the NE — SW direction is 
of order 0{N'^). However in most situations we believe that one can organize 
a set in the NE — SW direction substantially faster by applying the algorithm 
we shall now describe. 

Let us assume that the points {x{i),y{i)), i = 0, l,2,...,n have been orga- 
nized in the NE — SW direction. We shall now show how to find the parents 
of the next point {x{n + 1), y{n + 1)). 

Let mi be defined as the largest index for which 

x{m\) = max{a;(i) : x{i) < x{n + 1), y{i) < y{n + 1), < i < n} 
and let m2 be defined as the largest index for which 

t/(m2) = max{y(j) : x{j) < x(n + 1), y{j) <y{n+l), 0<i< n}. 

We now have two possibilities. Either mi = m2 or rui ^ m2- In the first 
case we are ready and (x(mi), y{m\)) is the only parent of {x{n + 1), y{n + 1)). 

If instead mi / m2 then we have obtained two parents to 
(.r(n-l-l), y{n+l)) namely (a;(mi), y{mi)) and (a;(m2), j/(m2)), and it is possible 
that there are further parents. In order to find these, we shall look for points 
{x{k),y{k)) e K satisfying 

a;(m2) < x{k) < x{mi) 

and 

y{mi) < y{k) < y{m2). 
We now define ms as the largest index for which 

x{m3) = max{a;(i) : a;(m2) < x{i) < x{mi), y{mi) < y{i) < y{m2), 1 < i < n}, 

if any such index exists, and we define 1714 as the largest index for which 

y{m4) = max{2/(j) : a;(m2) < x{j) < x{mi), y{mi) < y{j) < y{ni2), 1 < j < n}. 
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If ma is not defined, then there is no further parent. In case ms = rui then 
there is exactly one more parent namely (a;(m3), 2/(7713)). 

Otherwise we have found two new parents namely (^(tms), 7/(7713)) and (a;(7774), 7/(7714)). 
In this case there may be further parents namely if there are elements {x{k),y{k)) G 
K which satisfy 

x{m4) < x{k) < x{m3) 

and 

yirna) < y{k) < y{mi). 

To investigate whether there exist any points in this rectangle we proceed in 
the same way as above. 

We have now briefly described an algorithm by which one can organize the 
elements of Ki , K2 and Ki U K2 in such a way that it is easy to find the nearest 
elements in the NE — SW directions. 

Remark. As pointed out above it is easy to construct an example of a set with 
N pixels which has N'^/A parents and children. Whether this also requires that 
we need a storage for the structure of order N'^ is not 100 % sure, since in the 
more extreme situations it is likely that many "parent sets" are very similar, 
and therefore it is conceivable that these sets could be coded cflSciently. In the 
general random case we believe that the storage needed to store the information 
about parents and children will be much less than when N is large. 



14. SUMMARY. 

In this paper we have proved some results, which give rise to stopping criteria, 
when searching for new admissible arcs, when using the primal-dual algorithm 
for solving the transportation problem in the plane, in case the underlying 
distance-function d{-, •, •, •) is defined either by d{i,j,x,y) = {x — i)"^ + {y — j)^ 
or by 5{i,j, x, y) = \/{x — iY + {v — jY- We believe, that by using Theorems 5 
and 6 of Section 10, and Theorem 7 of section 12, it will be possible to reduce 
the computational complexity of the Euclidean transportation problem and the 
Euclidean assignment problem to approximately 0{N'^) from approximately 
0(iV^-^), which is the best limit so far, as far as the author knows. (See e.g. [2] 
and [14]). 
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